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Abstract 

This  paper  reports  on  the  computational  performance  and  detailed  results  of  ultra  large-scale  simulations  of  a  200  cm2  polymer  electrolyte 
fuel  cell  (PEFC)  using  a  23.5  million  gridpoint  mesh.  The  computer  code  is  based  on  a  comprehensive  single -phase  PEFC  model  that 
features  a  detailed  membrane-electrode  assembly  (MEA)  model,  electron  transport,  thermal  and  species  transport,  coolant  heat  transfer,  in 
addition  to  other  standard  functionalities.  Two  cases  under  dry  operation  are  simulated  and  compared.  One  case  concerns  an  infinitely  large 
coolant  flowrate  and  consequently  a  constant  temperature  of  bipolar  plates.  The  other  case  involves  a  finite  flowrate  and  a  lower  inlet  coolant 
temperature  designed  to  avoid  membrane  dryout  in  the  inlet  region  while  alleviating  electrode  flooding  in  the  outlet  region. 
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1.  Introduction 

Computer  simulation  of  PEFCs  has  been  an  active  research 
area:  see  the  most  recent  comprehensive  review  of  Wang  [1]. 
Significant  capabilities  exist  using  the  single-phase  model¬ 
ing  approach.  Notable  examples  include  a  detailed  model  of 
membrane-electrode  assembly  (MEA)  [2-4],  electron  trans¬ 
port  and  direct  prescription  of  constant  current  boundary 
condition  [5,6],  non-isothermal  modeling  and  hence  cou¬ 
pled  consideration  of  water  and  thermal  management  [7], 
variable  gas  flow  [8,9],  transient  analysis  [10],  and  experi¬ 
mental  validation  against  detailed  current  distribution  data 
[11,12]. 

For  single-phase  calculations,  Meng  and  Wang  [13]  first 
proposed  a  parallel  computing  methodology  to  handle  large- 
scale  simulations  involving  millions  of  gridpoints.  They  sim¬ 
ulated  a  variety  of  wet-to-dry  operating  conditions  with  a 
five-channel  serpentine  flow-field  and  benchmarked  the  par¬ 
allel  computing  performance  to  be  greater  than  7  x  speed-up 
with  10  processors.  The  largest  calculations  reported  so  far 
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were  by  Wang  [1]  (2.7  million  cells);  Shimpalee  et  al.  [14]  (5 
million  cells),  and  Wang  and  Wang  [15]  (2.7  million  cells). 
However,  the  industry  has  a  need  for  simulations  using  an 
order-of-magnitude  larger  mesh  with  reasonable  computing 
time.  This  paper  presents  23.5  million  gridpoint  calculations 
for  the  first  time,  herein  termed  ultra  large-scale  simulation  to 
distinguish  from  prior  studies  of  large-scale  simulations.  In 
addition,  we  will  show  a  capability  to  perform  co- simulation 
of  electrochemical/transport  phenomena  in  a  PEFC  and  heat 
transfer  in  coolant  channels.  Such  an  integrated  heat  and 
water  management  tool  is  highly  desirable  in  design  and  engi¬ 
neering  of  commercial-size  PEFCs. 


2.  Numerical  model 

Fig.  1  shows  the  computational  domain  of  the  200  cm2 
PEFC  featuring  a  24-channel,  3-pass  serpentine  flowfield, 
which  includes  the  flow  plates,  gas  and  coolant  channels, 
gas  diffusion  layers  (GDFs),  and  catalyst  layers  on  both  the 
anode  and  cathode,  as  well  as  the  membrane.  The  geometrical 
dimensions  of  various  components  in  the  cell  are  listed  in 
Table  1. 
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Nomenclature 

cp 

specific  heat  (Jkg-1  K-1) 

ck 

molar  concentration  of  species  k  (mol  m-3) 

D 

mass  diffusivity  of  species  (m2  s-1) 

F 

Faraday’s  constant  (96,487  C  eq._1) 

I 

current  density  (A  cm-2) 

K 

permeability  (m2) 

P 

pressure  (Pa) 

R 

gas  constant  (8 . 1 34  J  mol- 1  K- 1 ) 

S 

source  term  in  transport  equations 

T 

temperature  (K) 

u 

velocity  vector  (ms-1) 

Greek  letters 

a 

net  water  transport  coefficient  per  proton 

£ 

porosity 

<S> 

phase  potential  (V) 

K 

ionic  conductivity  (S  m-1) 

X 

membrane  water  content,  #H20/#S03- 

p 

density  (kgm-3) 

a 

electronic  conductivity  (S  cm-1) 

X 

shear  stress  (N  m-2) 

$ 

stoichiometric  flow  ratio 

Superscripts  and  subscripts 

a 

anode 

avg 

average 

c 

cathode 

e 

electrolyte 

eff 

effective  value 

g 

gas  phase 

k 

species 

m 

membrane 

s 

electronic  phase 

sat 

saturate  value 

w 

water 

2.1.  Governing  equations 


Thru- 

A 


plane 


Alony-channel  direction 


3-pass,  24-channel  Flowfield 


Coolant 

channel 


Bipolar  plate  Gas  channel  and  GDL  MEA 


Fig.  1.  Computational  domain  and  mesh  of  a  200  cm2  PEFC. 


(3) 

Energy  conservation  :  V  •  ( pcpuT )  =  V  •  (keiiVT)  +  Sj- 

(4) 

Proton  transport  equation  :  V  •  (/ce  v</>e)  +  S,pc  =  0 

(5) 

rr _ 

Electron  transport  equation  :  V  •  (cr  V0S)  +  5^=0 

(6) 

where  p,  u,  p,  Q,  T,  0e>  and  0S  denote  the  gas  density, 
superficial  fluid  velocity  vector,  pressure,  molar  concentra¬ 
tion  of  species  k ,  temperature,  electrolyte  and  electronic 
phase  potentials,  respectively.  The  species  considered  here 
are  hydrogen,  oxygen,  and  water.  The  various  source  terms, 
electrochemical  properties,  and  thermophysical  properties 
identified  for  various  regions  of  a  PEFC,  as  well  as  nec¬ 
essary  boundary  conditions,  have  been  detailed  in  [3-7] 
and  thus  are  not  repeated  here.  Select  physical  properties 
most  relevant  to  the  present  simulation  results  are  listed  in 
Table  1. 


The  single-phase  PEFC  model  employed  in  this  work 
consists  of  nonlinear,  coupled  partial  differential  equations 
governing  the  conservation  of  mass,  momentum,  species, 
energy,  and  charge  with  electrochemical  reactions.  The  equa¬ 
tions  can  be  written  in  the  vector  form,  as  [3-7]: 

Continuity  equation  :  V  •  ( pu )  =  0  (1) 

Momentum  conservation  : 

1 

-rV  •  ( puu )  =  —V p  +  V  •  r  +  Su  (2) 

£z 

Species  conservation  :  V  •  ( uCk )  =  V  •  {Df^VCp)  +  Sg 


2.2.  Numerical  procedure 

The  conservation  equations,  Eqs.  (1)— (6),  are  solved  by 
the  commercial  package,  Fluent®,  with  the  SIMPLEC  algo¬ 
rithm  [16],  using  a  parallel  computational  methodology  for 
a  Linux  PC  cluster.  The  SIMPLEC  algorithm  is  to  update 
the  pressure  and  velocity  fields  from  the  solution  of  the 
pressure  corrections  equation,  solved  by  algebraic  multi¬ 
grid  (AMG)  method  [17].  Following  the  solution  of  the  flow 
field,  the  energy,  species,  proton,  and  electron  equations  are 
solved.  Approximately  23.5  million  computational  elements 
(46  x  900  x  600)  are  employed  to  capture  details  of  flow, 
thermal,  electrochemical  and  mass  transport  phenomena  in 
the  PEFC.  The  MEA,  where  the  important  electrochemical 


132 


Y.  Wang,  C.-Y.  Wang  /  Journal  of  Power  Sources  153  (2006)  130-135 


Table  1 


Geometrical  and  physical  parameters  [9,19] 


Quantity 

Value 

Gas  channel  depth 

1.0  mm 

Gas  channel  width 

1.0  mm 

Land  width 

1.0  mm 

Diffusion  layer  thickness 

0.2  mm 

Catalyst  layer  thickness 

0.01  mm 

Membrane  thickness 

0.045  mm 

Fuel  cell  footprint  area 

143  mm  x  143  mm 

Anode/cathode/coolant  flow 

Co-flow 

Anode/cathode  pressure,  P 

2. 0/2.0  atm 

Average  current  density,  7av g 

1.0  Acm~2 

Stoichiometric  ratio  £  in  the  anode/cathode 

2. 0/2.0 

RH  of  anode/cathode  inlet  @  80  °C 

50%/50% 

GDL  porosity,  e 

0.6 

Catalyst  layer  porosity,  sg 

0.4 

Volume  fraction  of  ionomer  in  catalyst  layer,  sm 

0.26 

GDL  permeability,  K 

10“12m2 

Specific  heat,  cp,  of  coolant  (de-ionized  water) 

4182Jkg-1  K-1 

Water  activity,  a 

a  =  Cwsf,T ,  where  logio  PSM  =  2. 1794  +  0.02953  (T  273. 15) 

-  9.1837  x  10"5(r-  273.15)2  +  1.4454  x  10"7(r-  273.15)3 

Water  content  in  membrane,  k 

J  0.043  +  17.81a  —  39.85a2  +  36.0a3  for  0  <  a  <  1 
=  1  14+  1.4(a-  1)  for  1  <a  <  3 

reactions  and  water  back-diffusion  process  occur,  is  accu¬ 
rately  captured  by  placing  five  and  eight  gridpoints  in  the 
catalyst  layer  and  the  membrane,  respectively,  along  the 
through-plane  direction.  Details  of  the  computational  mesh 
are  shown  in  Fig.  1.  A  constant  average  current  density  of 
1  A  cm-2  is  specified  as  an  input  parameter,  allowing  the  local 
current  density  and  electronic  phase  potential  to  be  predicted 
to  vary  spatially  according  to  local  operating  conditions  [6]. 
A  simulation  using  23.5  million  cells  typically  requires  about 
600  iterations  with  the  global  mass  balance  less  than  1%  and 
residuals  of  species  equations  less  than  10-6,  and  takes  nearly 
20  h  on  32  nodes  of  2.8  GHz  Intel  Pentium  4  CPU  and  1 .0  GB 
DDR  RAM. 

3.  Results  and  discussion 

A  200  cm2  PEFC  is  considered  with  the  inlet  gases  having 
relative  humidity  (RH)  of  50%  relative  to  80  °C  on  both  anode 
and  cathode  sides.  Other  common  operating  conditions  are 
detailed  in  Table  1.  Two  simulation  cases  are  presented  in 
this  work.  Case  1  involves  an  infinitely  large  flow  rate  of  the 
coolant,  thus  effectively  keeping  the  flow  plate  temperature 
constant  at  80  °C  or  353.15  K.  In  contrast,  Case  2  simulates  a 
thermal  management  scheme  to  enhance  cell  performance  in 
which  the  coolant  temperature  is  selected  at  72  °C  at  the  inlet 
and  its  flowrate  at  2.9  x  10-6  ms-1.  The  lower  inlet  coolant 
temperature  helps  to  increase  the  local  gas  RH  near  the  inlet 
and  hence  better  hydrates  the  membrane,  whereas  the  small 
coolant  flowrate  is  used  to  allow  the  coolant  temperature  to 
increase  by  ~10  °C  from  the  inlet  to  outlet  so  as  to  alleviate 
possible  cathode  flooding  in  the  exit  region,  as  well  as  to  ele¬ 
vate  the  coolant  exit  temperature  for  better  heat  dissipation  in 


the  radiator.  The  resulting  cell  voltages  predicted  under  the 
two  cooling  conditions  are  contrasted  in  Table  2.  There  is  a 
26  mV  gain  in  cell  voltage  in  Case  2  as  compared  with  Case 
1 .  In  addition,  detailed  contours  of  the  current  density,  mem¬ 
brane  temperature,  membrane  water  content,  and  membrane 
water  crossover  coefficient  are  presented  below  to  illustrate 
the  predictive  capabilities  of  such  an  ultra  large-scale  simu¬ 
lation  tool. 

Fig.  2  shows  the  current  density  distribution  in  Case  1.  It 
can  be  seen  that  the  current  density  initially  increases  along 
the  cathode  flow,  as  can  be  explained  by  the  ohmic  control  of 
cell  performance.  The  dry  anode  and  cathode  gases  result  in 
dryout  of  the  membrane  in  the  inlet  region,  thereby  increas¬ 
ing  the  local  ionic  resistance.  As  gases  travel  through  the 
fuel  cell,  membrane  hydration  improves  by  product  water 
and  hence  its  ionic  resistance  decreases  downstream  and  the 
current  density  increases. 

Fig.  3  shows  the  temperature  distribution  in  the  mem¬ 
brane  in  Case  1.  Despite  that  the  current  collecting  lands 
maintain  the  constant  temperature  of  80  °C  due  to  infinitely 
large  coolant  flowrate,  there  is  ~3  K  temperature  rise  in  the 
membrane  due  to  the  limited  in-plane  thermal  resistance  of 
GDL.  In  addition,  the  membrane  temperature  increases  along 
the  cathode  flow,  in  accordance  with  the  increasing  current 
density. 

Table  2 


Coolant  conditions  and  cell  voltages 


Case  1 

Case  2 

Coolant  flow  rate 

Infinite 

2.9  x  10-6 ms-1 

Coolant  inlet  temperature 

353.15  K  (80  °C) 

345.15  K  (72  °C) 

Predicted  cell  voltage 

0.590  V 

0.616V 
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Fig.  2.  Current  density  distribution  in  the  membrane  in  Case  1. 


Fig.  4.  Water  content  at  the  membrane-cathode  interface  in  Case  1. 


Fig.  4  presents  the  distribution  of  membrane  water  content 
at  the  interface  between  the  membrane  and  cathode  catalyst 
layer  for  Case  1 .  It  can  be  seen  that  the  cathode  is  dry  in  the 
inlet  region  in  spite  of  water  production.  In  addition,  water 
content  reaches  the  fully  hydrated  level,  i.e.  A  =  14,  in  the  out¬ 
let  region,  resulting  in  likely  occurrence  of  cathode  flooding. 

Fig.  4  points  out  two  problems  inherent  in  the  design  of 
Case  1 .  First,  the  dry  gas  feed,  while  preferred  for  the  automo¬ 
tive  application,  causes  membrane  dryout  in  the  inlet  region. 
Secondly,  the  outlet  gas  stream  is  over-saturated.  These  two 
issues  can  be  resolved  by  proper  thermal  control  through  the 
coolant  flow  design,  such  as  in  Case  2.  Fig.  5  presents  the  tem¬ 
perature  distribution  in  the  cathode  flow  plate  for  Case  2.  It  is 
seen  that  this  coolant  design,  using  a  lower  inlet  temperature 
and  a  small  flowrate,  creates  a  thermal  gradient  (345-355  K) 
from  the  inlet  to  the  outlet  regions.  A  significant  benefit  of 
this  thermal  gradient  along  the  flow  direction  can  be  seen  in 
Fig.  6,  where  membrane  water  content  is  displayed.  As  com¬ 
pared  to  Fig.  3  with  isothermal  cooling,  Case  2  increases  the 


membrane  water  content  in  the  inlet  region  under  lower  local 
temperatures  and  at  the  same  time  maintains  fully  saturated 
gas  in  the  outlet  region  under  higher  temperatures.  The  net 
result  of  this  thermally  graded  fuel  cell  is  that  the  membrane 
water  content  becomes  more  uniform. 

As  a  result  of  the  more  uniform  degree  of  membrane 
hydration,  the  current  density  near  the  inlet  becomes  higher, 
as  shown  in  Fig.  7.  In  addition,  more  uniform  current  distri¬ 
bution  results. 

It  is  worth  noting  that  the  coolant  condition  in  Case  2,  how¬ 
ever,  results  in  greater  temperature  gradient  in  the  membrane, 
as  shown  in  Fig.  8.  There  is  ~12K  temperature  difference 
between  the  inlet  and  outlet  areas.  This  temperature  non¬ 
uniformity  could  potentially  cause  membrane  degradation 
due  to  thermal  stresses  [18]. 

Fig.  9  displays  the  net  water  transfer  coefficient  through 
the  membrane,  a ,  for  Case  2.  It  can  be  seen  that  this 
coefficient,  indicating  water  crossover  through  the  mem¬ 
brane,  varies  spatially.  In  the  inlet  region  of  both  anode  and 
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Fig.  3.  Membrane  temperature  distribution  in  Case  1. 


Fig.  5.  Surface  temperature  at  the  cathode  flow  plate  in  Case  2. 
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Fig.  6.  Water  content  at  the  membrane-cathode  interface  in  Case  2. 


Fig.  9.  Contours  of  net  water  transfer  coefficient,  a,  across  the  membrane 
in  Case  2. 


Fig.  7.  Current  density  distribution  in  the  membrane  in  Case  2. 
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Fig.  8.  Membrane  temperature  contours  in  Case  2. 


cathode,  the  anode  loses  water  to  the  cathode.  Once  the  cath¬ 
ode  gas  gains  moisture  from  product  water,  the  net  water  flow 
is  directed  from  the  cathode  to  anode  due  to  dominance  of 
water  back-diffusion. 


4.  Conclusions 

Ultra  large-scale  simulations  of  a  commercial- size  PEFC 
with  24-channel  3 -pass  flowfield  have  been  carried  out  to 
demonstrate  the  computational  performance  of  a  parallelized 
code  for  PEFCs,  and  to  illustrate  the  internal  phenomena  and 
distributions  of  key  parameters  in  a  representative  industrial 
cell  as  well  as  the  profound  coolant  effects.  Using  a  32- 
processor  PC  cluster,  we  have  shown  a  computing  time  of 
roughly  20  h  for  a  simulation  using  a  23.5  million  gridpoint 
mesh  that  considers  comprehensive  single-phase  phenomena 
in  a  PEFC.  In  addition,  the  computed  contours  of  current 
and  water  distributions  indicate  that  low-humidity  operation 
commonly  used  in  automotive  applications  leads  to  mem¬ 
brane  dryout  in  the  inlet  region  and  cathode  flooding  in 
the  outlet  region.  These  two  seemingly  conflicting  issues 
can  be  simultaneously  resolved  by  using  a  cooling  scheme 
where  a  lower  inlet  coolant  temperature  and  a  small  coolant 
flowrate  are  employed  to  create  a  thermal  gradient  along  the 
cathode  gas  flow.  The  resulting  “cold”  inlet  region  main¬ 
tains  high  relative  humidity  despite  the  dry  gas  supply,  and 
the  “hot”  outlet  region  alleviates  susceptibility  of  cathode 
flooding. 

While  the  present  paper  exemplifies  only  one  application 
of  the  ultra  large-scale  simulation  tool,  numerous  other  appli¬ 
cations  of  the  same  computational  intensity  have  been  carried 
out  routinely  at  this  research  laboratory  as  well  as  several 
automotive  companies. 
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